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Abstract 


An important problem that arises in different areas of science and engineering is 
that of computing the limits of sequences of vectors {x m }, where x m £ C N , N 
being very large. Such sequences arise, for example, in the solution of systems of 
linear or nonlinear equations by fixed-point iterative methods, and lim, Jl _ ) . 00 x m 
are simply the required solutions. In most cases of interest, however, these se¬ 
quences converge to their limits extremely slowly. One practical way to make the 
sequences { x m } converge more quickly is to apply to them vector extrapolation 
methods. Two types of methods exist in the literature: polynomial type methods 
and epsilon algorithms. In most applications, the polynomial type methods have 
proved to be superior convergence accelerators. Three polynomial type meth¬ 
ods are known, and these are the minimal polynomial extrapolation (MPE), the 
reduced rank extrapolation (RRE), and the modified minimal polynomial extrap¬ 
olation (MMPE). In this work, we develop yet another polynomial type method, 
which is based on the singular value decomposition, as well as the ideas that lead 
to MPE. We denote this new method by SVD-MPE. We also design a numerically 
stable algorithm for its implementation, whose computational cost and storage re¬ 
quirements are minimal. Finally, we illustrate the use of SVD-MPE with numerical 
examples. 


Mathematics Subject Classification 2000: 15A18, 65B05, 65F10, 65F50, 65H10. 


Keywords and expressions: Vector extrapolation, minimal polynomial extrapola¬ 
tion, singular value decomposition, Krylov subspace methods. 



1 Introduction and background 


An important problem that arises in different areas of science and engineering is that 
of computing limits of sequences of vectors {* m }0 where x m € C^, the dimension 
N being very large in many applications. Such vector sequences arise, for example, 
in the numerical solution of very large systems of linear or nonlinear equations by 
fixed-point iterative methods, and lim m _^ 00 x m are simply the required solutions to 
these systems. One common source of such systems is the finite-difference or finite- 
element discretization of continuum problems. In later chapters, we will discuss further 
problems that give rise to vector sequences whose limits are needed. 

In most cases of interest, however, the sequences {x m } converge to their limits ex¬ 
tremely slowly. That is, to approximate s = linim^oo x m , with a reasonable prescribed 
level of accuracy, by x m , we need to consider very large values of m. Since, the vectors 
x m are normally computed in the order m = 0,1, 2,... , it is clear that we have to 
compute many such vectors until we reach one that has acceptable accuracy. Thus, 
this way of approximating s via the x m becomes very expensive computationally. 

Nevertheless, we may ask whether we can do something with those x m that are 
already available, to somehow obtain new approximations to s that are better than 
each individual available x m . The answer to this question is in the affirmative for 
at least a large class of sequences that arise from fixed-point iteration of linear and 
nonlinear systems of equations. One practical way of achieving this is by applying 
to the sequence {x m } a suitable convergence acceleration method (or extrapolation 
method). 

Of course, in case linim^oo x m does not exist, it seems that no use could be made of 
the x m . Now, if the sequence {x m } is generated by an iterative solution of a linear or 
nonlinear system of equations, it can be thought of as “diverging from” the solution s 
of this system. We call s the antilimit of {x m } in such a case. It turns out that vector 
extrapolation methods can be applied to such divergent sequences {x m } to obtain good 
approximations to the relevant antilimits, at least in some cases. 

Two different types of vector extrapolation methods exist in the literature: 

1. Polynomial type methods: The minimal polynomial extrapolation (MPE) of Cabay 
and Jackson [3], the reduced rank extrapolation (RRE) of Eddy [5] and Mesina 
m, and the modified minimal polynomial extrapolation (MMPE) of Brezinski 
[3], Pugachev [13] and Sidi, Ford, and Smith [23] . 

2. Epsilon algorithms: The scalar epsilon algorithm (SEA) of Wynn [28] (which 
is actually a recursive procedure for implementing the transformation of Shanks 
m ), the vector epsilon algorithm (VEA) of Wynn [29], and the topological epsilon 
algorithm (TEA) of Brezinski [3]. 

The paper by Smith, Ford, and Sidi [25] gives a review of all these methods (ex¬ 
cept MMPE) that covers the developments in vector extrapolation methods until the 
end of the 1970s. For up-to-date reviews of MPE and RRE, see Sidi [20] and [21]. 

1 Unless otherwise stated, {c m } will mean {c m }m=o throughout this work. 


1 



Numerically stable algorithms for implementing MPE and RRE are given in Sidi [ T8] , 
these algorithms being also economical both computationally and storagewise. For the 
convergence properties and error analyses of MPE, RRE, MMPE, and TEA, as these 
are applied to vector sequences generated by fixed-point iterative methods from linear 
systems, see the works by Sidi m , m, [19] . Sidi, Ford, and Smith [23], Sidi and 
Bridger [22], and Sidi and Shapira [23]. VEA has been studied by Brezinski ra> 0 , 
Gekeler [6], Wynn [30] , m, and Graves-Morris m, uni. 

Vector extrapolation methods are used effectively in various branches of science 
and engineering in accelerating the convergence of iterative methods that result from 
large sparse systems of equations. 

All of these methods have the useful feature that their only input is the vector 
sequence {x m } whose convergence is to be accelerated, nothing else being needed. 
In most applications, however, the polynomial type methods, especially MPE and 
RRE, have proved to be superior convergence accelerators; they require much less 
computation than, and half as much storage as, the epsilon algorithms for the same 
accuracy. 

In this work, we develop yet another polynomial type method, which is based on the 
singular value decomposition (SVD), as well as some ideas that lead to MPE. We denote 
this new method by SVD-MPE. We also design a numerically stable algorithm for its 
implementation, whose computational cost and storage requirements are minimal. The 
new method is described in the next section. In Section [3] we show how the error in 
the approximation produced by SVD-MPE can be estimated at zero cost in terms of 
the quantities already used in the construction of the approximation. In Section [4] 
we give a very efficient algorithm for implementing SVD-MPE. In Section [5] we derive 
determinant representations for the approximations produced by SVD-MPE, while in 
Section [6] we show that this method is a Krylov subspace method when applied to 
vector sequences that result from the solution of linear systems via fixed-point iterative 
schemes. Finally, in Section [7] we illustrate its use with two numerical examples. 

Before closing, we state the (reduced version of) the well known singular value 
decomposition (SVD) theorem. For different proofs, we refer the reader to Golub and 
Van Loan [8], Horn and Johnson HD. Stoer and Bulirsch [26], and Trefethen and Bau 
[27J, for example. 

Theorem 1.1 Let A € C rxs , r > s. Then there exist unitary matrices U € C rxs , 
V E C sxs , and a diagonal matrix £ = diag(ai,... ,<r s ) E M sxs , with o\ > cr 2 > ■ ■ ■ > 
a s > 0, such that 

A = U'EV*. 

Furthermore, if U = [ u\ \ ■ ■ ■ \ u s ] and V = [rq | • • • | v s ], then 

Avi = OiUi , A*u.i = OiVi , i = 1,..., s. 

In case rank( A) = t, there holds <Tj > 0, i = 1,..., t, and the rest of the are zero. 

Remark: The ai are called the singular values of A and the and u, are called the 

corresponding right and left singular vectors of A, respectively. We also have 

A*Avi = afvi, AA*Ui = ofiij, i = 1,..., s. 
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2 Development of SVD-MPE 


In what follows, we use boldface lower case letters for vectors and boldface upper case 
letters for matrices. In addition, we will be working with general inner products (•, •) 
and the I 2 norms || • || induced by them: These are defined as follows: 

• In C N , with M E C NxN hermitian positive definite, 

(a, b) M = a*Mb, \\a\\ M = V( a ,a) M • (2.1) 

• In C fc+1 , k = 1,2,..., with L k E ([4 fc+1 i x ( fc+1 ) hermitian positive definite, 



Of course, the standard Euclidean inner product a*b and the I 2 norm \/a*a induced 
by it are obtained by letting M = I in (12.1)1 and L k = I in (12.21) : we will denote these 
norms by || • || 2 - (We will denote by I the identity matrix in every dimension.) 

2.1 Summary of MPE 

We begin with a brief summary of minimal polynomial extrapolation (MPE). We use 
the ideas that follow to develop our new method. 

Given the vector sequence {x m } in C^, we define 

— ®m+l ^ — 0,1,..., (2-3) 

and, for some fixed n, define the matrices U k via 

U k = [u n | u n+1 | ••• | u n+k } € C N ^ k+1 \ (2.4) 

Clearly, there is an integer ko < N, such that the matrices U k , k = 0,1,..., ko — 1, 
are of full rank, but U ko is not; that is, 

rank (U k ) = k + 1, k = 0,1,..., ko — 1; rank (U k 0 ) = ko. (2-5) 

(Of course, this is the same as saying that {wo, «i, ■ ■ ■, u ko -i} is a linearly independent 
set, but {wo, u ii ■ ■ ■ i u k 0 } is not.) Following this, we pick a positive integer k < ko and 
let the vector d = [cq, ci, ..., c k - i] T S be the solution to 


min 

Co,Cl,...,Cfc_i 


fc-i 


^ ( CjU. n _|_j T u n j rk 


i=0 


M 


This minimization problem can also be expressed as in 


( 2 . 6 ) 


min \\U k -ic' + u n+k \\ M , (2.7) 

C 
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and, as is easily seen, c' is the standard least-squares solution to the linear system 
Uk-ic' = —u n+ k, which, when k < N, is overdetermined, and generally inconsistent. 
With co, ci,..., Cfc_i determined, set c k = 1, and compute the scalars 70 , 71 ,..., 7 k via 

7 i = —^-, i = 0,l,...,k, ( 2 . 8 ) 

£ i= o c i 

provided Xj=o c j 7 ^ 0 - Note that 


k 


= L 


(2.9) 


Finally, set 


k 

Sn,k = ^ ^ Ti®n-|-i 
i =0 


( 2 . 10 ) 


as the approximation to s, whether s is the limit or antilimit of {x m }. 

What we have so far is only the definition (or the theoretical development) of MPE 
as a method. It should not be taken as an efficient computational procedure (algo¬ 
rithm), however. For this topic, see {18] , where numerically stable and computationally 
and storagewise economical algorithms for MPE and RRE are designed for the case in 
which M = I. A well documented FORTRAN 77 code for implementing MPE and 
RRE in a unified manner is also provided in {T 8 , Appendix B]. 


2.2 Development of SVD-MPE 

We start by observing that the unconstrained minimization problem for MPE given in 
m can also be expressed as a superficially “constrained” minimization problem as 
in 

min \\U k c\\ M , subject to c k = 1; c = [c 0 , ci,..., c k ] T . (2.11) 

C 

For the SVD-MPE method, we replace this “constrained” minimization problem by 
the following actual constrained minimization problem: 

min \\Ukc\\ M , subject to \\c\\ L =1- c = [c 0 , ci,..., c k ] T ■ (2.12) 

C 

With Co, ci,..., c k determined, we again compute 7 oj 71 , • • •, 7fc via 

7i = —^-, i = 0,1,..., k, (2.13) 

X;=o c i 

provided Ylj =0 c j 7^ 0, noting again that 

k 

X> = L ( 2 - 14 ) 

i=0 
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Finally, we set 


k 

&n,k ^ ' T i^n+i (2.15) 

i =0 

as the SVD-MPE approximation to s, whether s is the limit or antilimit of {x m }. 

Of course, the minimization problem in (12.121) has a solution for c = [co, ci,..., Cfc] T . 
Let f7 m in = \\U k c\\ M for this (optimal) c. Lemma lTTI that follows next gives a complete 
characterization of cr m [ n and the (optimal) c. 


Lemma 2.1 Let (T^o, a k i, ■ ■ ■, crkk be the singular values of the N x 

(k + 1) matrix 

U k = MV 2 U k L~ 1/2 , 

(2.16) 

ordered as in 

O'kO > CTfcl > ■ ■ ■ > &kk, 

(2.17) 

and let h k i be the corresponding right singular vectors ofU k , that is, 

UfcUfahki &ki^ J ki’) ||^fci||2 1? ^ 0,1, . . . , /c. 

(2.18) 

Assuming that a k k, the smallest singular value ofU k , is simple, the (optimal) solution 
c to the minimization vroblem in (12.120 is uniaue fuv to a multivlicative constant r. 
r =1), and is given as in 

C — Lfc ^min — — &kk' 

(2.19) 

~ 1/2 

Proof. The proof is achieved by observing that, with c = Lf c, 


\\U k c\\ M = \\U k c \\2 and \\c\\ Lk = ||c|| 2 , 

(2.20) 

so that the problem in (12.121) becomes 


nun \\U k c 2 , subject to |c 2 = 1- 

C 

(2.21) 

We leave the details to the reader. 

■ 


Remarks: 

1. In view of the nature of the solution for the (optimal) c involving singular values 
and vectors, as described in Lemma 12.11 we call this new method SVD-MPE. 

2. Note that if rank (Uk) = k + 1, then rank ([/*,) = k + 1 too, and we therefore 
have akk > 0 . 

3. Of course, exists if and only if the (optimal) c = [co, c±,... , Cfc] T satisfies 
^2j=q Cj / 0. In addition, by (12.131) . the 7 j are unique when a k k is simple. 

Before we go on to the development of our algorithm in the next section, we state 
the following result concerning the finite termination property of SVD-MPE, whose 
proof is very similar to that pertaining to MPE and RRE given in ED: 
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Theorem 2.2 Let s be the solution to the nonsingular linear system x = Tx-\-d, and 
let {x m } be the sequence obtained via the fixed-point iterative scheme x m+ i = Tx m +d, 
m = 0,1,, with xq chosen arbitrarily. If k is the degree of the minimal polynomial 
of T with respect to e n = x n — s (equivalently, with respect to w n jjl then s n:k produced 
by SVD-MPE satisfies s n>k = s. 


3 Error estimation 

We now turn to the problem of estimating at zero cost the error s n ^ — s, whether s is 
the limit or antilimit of {x m }. Here we assume that s is the solution to the system of 
equations 

* = /(*); f :C N ^C N , xe C N , 

and that the vector sequence {x m } is obtained via the fixed-point iterative scheme 

x m -\-1 — f{p^rn)i tTl — 0 , 1 ,..., 

Xq being the initial approximation to the solution s. 

Now, if x is some approximation to s, then a good measure of the error x — s in 
x is the residual vector r{x) of x, namely, 

r(x) = f(x) - x. 

This is justified since lima,-^ r(x) = r(s) = 0. We consider two cases: 

i- f (*) is linear; that is, f{x) = Tx+d, where T G (£ NxN anc [ I—T is nonsingular. 
In this case, we have 

r(x) = Tx + d — x = (T — /)(* — s), 

and, therefore, by Ya=qH = s n,k = Y)i=oH x n+i satisfies 

k k 

= ^ ^ 'fiiiTXn+i d) ^n+z] = ^ 1 ® n+i ) 

2=0 2=0 

and thus 

r(s nik ) = U k 7, 7 = [70,71,... ,7fe] T . 

2. f(x) is nonlinear. 

In this case, assuming that \im. in ^ f00 x rn = s and expanding f(x m ) about s, we 
have 

X m + 1 = f{s) + F(s)(x m - s) + 0(||£C m - s|| 2 ) as m oo, 

2 Given a matrix B € C rxr and a nonzero vector a G C r , the monic polynomial P( A) is said to be 
a minimal polynomial of B with respect to a if P(B)a = 0 and if P( A) has smallest degree. 

It is easy to show that the minimal polynomial P( A) of B with respect to a exists, is unique, and 
divides the minimal polynomial of B, which in turn divides the characteristic polynomial of B. [Thus, 
the degree of P( A) is at most r, and its zeros are some or all of the eigenvalues of _B.[ Moreover, if 
Q{B)a = 0 for some polynomial Q( A) with degQ > degP, then P( A) divides Q( A). Concerning this 
subject, see Householder m, for example. 


i =0 


''l ] r )i u n+ii 

(3.1) 
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where F(x) is the Jacobian matrix of the vector-valued function fix) evaluated 
at x. Recalling that s = f(s), we rewrite this in the form 

x m+1 = s + F(s){x m - s) + 0{\\x m - s|| 2 ) as m -)► oo, 

from which, we conclude that the vectors x m and x m+ i satisfy the approximate 
equality 

x m+ i « s + F{s){x m — s ) for all large m. 

That is, for all large rn. the sequence {x m } behaves as if it were being generated 
by an IV-dimensional approximate linear system of the form ( I—T)x ~ d through 

Xm+l ~ Fx m + d, TTl — 0 , 1 , , 

where T = F[s) and d = [I — F{s)]s. In view of what we already know about 
r i s n,k) for linear systems [from (13.11) ]. for nonlinear systems, close to convergence, 
we have 

r(s Htk ) 7 = bin 7ii ■ ■ • >7fc] T - (3-2) 

Now, we can compute ||£/*; 7 ||m at no cost in terms of the quantities that result 
from our algorithm, without having to actually compute U k 7 itself. Indeed, we have 
the following result: 

Theorem 3.1 Let a kk be the smallest singular value of U k and let h kk be the corre¬ 
sponding right singular vector. Then the vector U k 7 resulting from s U)k satisfies 

\\U kl \\ M = ,° kk (3.3) 

l£i=o c il 


Proof. First, the solution to ( 12 . 121 ) is c = L l ^ 2 h kk by (12.191) . Next, letting a = 
^y = 0 Cj, we have 7 = c/a by (12.131) . Consequently, 


Thus, by Lemma l2.ll we have 


U k 7 


UkC 

a 


||C7 fe7ll JW 


jjt/fegjjM 

lal 


@kk 

lal 


which is the required result. 


4 Algorithm for SVD-MPE 

We now turn to the design of a good algorithm for constructing numerically the ap¬ 
proximation s nk that results from SVD-MPE. We note that matrix computations in 
floating-point arithmetic must be done with care, and this is what we would like to 
achieve here. 
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In this section, we let M = I and L k = I for simplicity. Thus, U k = Uk- Since 
there is no room for confusion, we will also use <7j, hi, and g i to denote a k i, h k i, and 
9kii respectively. 

As we have seen in Section [2J to determine s nyk , we need h k , the right singular 
vector of Uk corresponding to its smallest singular value a k . Now, a k and h k can 
be obtained from the singular value decomposition (SVD) of U k € (^Nx(k+ 1 ). of 
course, the SVD of U k can be computed by applying directly to Uk the algorithm of 
Golub and Kahan [Tj, for example. Here we choose to apply SVD to Uk in an indirect 
way, which will result in a very efficient algorithm for SVD-MPE that is economical 
both computationally and storagewise in an optimal way. Here are the details of the 
computation of the SVD of Uk, assuming that rank (U k ) = k + 1: 


1. We first compute the QR factorization of Uk in the form 

U k = Q k R k ; Q k g c Nx ( k+1 \ R k e c( fc+1 ) x ( fc+1 ), (4.1) 

where Q k is unitary (that is, Q* k Q k = I ) and R k is upper triangular with positive 
diagonal elements, that is, 


Qk = [Qo\Qi \ ■■■ I Qkl R k = 


r oo R)i 

m 


rok 

r\k 

Tkk 


q*qj = 5ij V i, j; r tj = q*Uj V i < j\ ru> 0 Vi. 


(4.2) 

(4.3) 


Of course, we can carry out the QR factorizations in different ways. Here we do 
this by the modified Gram-Schmidt process (MGS) as follows: 


1. Compute ?’oo = ||wn|| 2 and q 0 = u n /r 00 . 

2. For j = 1,. .., k do 

Set up = u n+ j 

For i = 0 , 1,..., j — 1 do 

* M i (i+i) M 

rij = q\ u\ and u\ = uy - r ij q i 

end do (i) 

Compute rjj = ||w ^||2 and qj = up/rjj. 

end do ( j) 


Note that the matrices Q k and R k are obtained from Q k _\ and R k -i, respec¬ 
tively, as follows: 





rok 

Qk — [Qk- 11 9k]i 

Rk = 

Rk— 1 

T’k—l,k 



i 

O 

o 

Tkk 


(4.4) 


For MGS, see [8], for example. 








2. We next compute the SVD of R k : By Theorem ll.il there exist unitary matrices 
Y,H £ c (fc+1)x(fe+1) , 

Y = [y 0 \y 1 \---\y k ], H=[h 0 \h 1 \ ■■■ \ h k ]; Y*Y = I , H*H = I 

(4.5) 

and a square diagonal matrix S € M( fc + 1 ) x ( fe + 1 ) ; 

S = diag (<r 0 , < 7 i,..., <Tfc); cr 0 > 07 > • • • > cr fc > 0, (4.6) 

such that 

R k = YYH*. (4.7) 

In addition, since R k is nonsingular by our assumption that rank (U k ) = k + 1, 
we have that a % > 0 for all i. Consequently, a k > 0. 

3. Substituting (14.71) in (14.11) . we obtain the following true singular value decompo¬ 
sition of U k : 

U k = GYH*- G = Q k Y € C N ^ k+l ) unitary, G*G = I; ^ 

G = [g 0 \gi I ••• \g k ], 9i9j = <%• 

Thus, ai, the singular values of R k , are also the singular values of U k , and 
hi , the corresponding right singular vectors of R k , are also the corresponding 
right singular vectors oiU k . [Of course, the are corresponding left singular 
vectors of U k . Note that, unlike Y, H , and X, which we must compute for our 
algorithm, we do not need to actually compute G. The mere knowledge that the 
SVD of U k is as given in (14.81) suffices to conclude that c = h k is the required 
optimal solution to ( 12 . 121 ) . and continue with the development of our algorithm.] 


With c = h k already determined, we next compute the 7 i as in (12.131) : that is, 

7 = ( 4 - 9 ) 

E j= o c j 

provided Ylj=o c j / °- 
Next, by the fact that 

i —1 

•Kn+i — T ^ ) ^n+j 

3= 0 

and by (12.141) . we can re-express s n ^ k in (12.151) in the form 
k -1 

Sn,k = T ^ ^ ^j^n+j = *n T U k —i£] £ = [£o> £l; • • ■ > l] j (4-10) 

3=0 

where the are computed from the 7 j as in 

k 

f-i = i; ij = ^2 n = £?-i -7 31 J = 0 , 1 , ■ ■ ■, k - 1- ( 4 - n ) 

i=j+l 
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Making use of the fact that U k -i = Q k _ k Rk-i, with 


Qk-1 = [Qo\Qi 


Qk-i. 


Rk—1 — 


r oo Rn 

rn 


r 0 ,k-i 

n,k-i 

T’k—l,k—l 


(4.12) 


where the q i and the rjj are exactly those that feature in (14. 2p and ()4.3p . we next 
rewrite (14.1011 as in 

Sn,k = 4“ Qk— l(-^fc—l£)- (4.13) 

Thus, the computation of s n ^ k can be carried out economically as in 

S n ,k = x n + Qk-iV, V = Rk- i£, V = [vo,Vi,-■ ■ ,Vk-i] T ■ (4.14) 


Of course, Q k _iT) is best computed as a linear combination of the columns of Q k _ 1; 
hence (|4.14p is computed as in 

k -1 

s n ,k = x n + ^2rnq i . (4.15) 

i=o 

It is clear that, for the computation in ([4.14ft and (14.151) . we need to save both Q k and 

Rk- 

This completes the design of our algorithm for implementing SVD-MPE. For con¬ 
venience, we provide a systematic description of this algorithm in Table 14.11 

Note that the input vectors x n+ i, i = l,...,k + 1, need not be saved; actually, 
they are overwritten by u n+ i, i = 0,1 ,k, as the latter are being computed. As 
is clear from the description of MGS given above, we can overwrite the matrix U k 
simultaneously with the computation of Q k and R k , the vector q n+3 overwriting u n+3 
as soon as it is computed, j = 0 , 1 ,..., k\ that is, at any stage of the QR factorization, 
we store k +2 A r -dimensional vectors in the memory. Since N » k in our applications, 
the storage requirement of the (k + 1) x (k + 1) matrix R k is negligible. So is the cost 
of computing the SVD of R k , and so is the cost of computing the (k + l)-dimensional 
vector r/. Thus, for all practical purposes, the computational and storage requirements 
of SVD-MPE are the same as those of MPE. 


Remark: If we were to compute the SVD of U k , namely, U k = GYiH'. directly, 
and not by (i) first carrying out the QR factorization of U k as Uk = Q k Rk, and 
(ii) then computing the SVD of R k as R k = VSiT*, then we would need to waste 
extra resources in carrying out the computation of s Utk = Ya=o li x n+i = x n + Uk- 1 £. 

1. We would either have to save U k while computing the matrix G in its singular 
value decomposition. Thus, we would need to save two N x (k + 1) matrices, 
namely, U k and G in core memory simultaneously. 

2. In case we are worried about storage, therefore, do not wish to save Uk, we need 
to recompute the vectors x n , cc n+ i,..., x n+k in order to compute s n k- Thus, we 
would need to compute these vectors twice to complete the determination of s n ^- 

Thus, the approach we have proposed here for carrying out the singular value decom¬ 
position of Uk enables us to save extra computing and storage very conveniently. 
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Step 0 . Input: The vectors x n , x n +i,..., x n +k+i- 

Step 1. Compute Ui = A Xi = Xi+ 1 — Xi, i = n, n + 1,..., n + k. 
Set U k = [u n | u n+ 1 1 • • • | u n+k ] e C Arx(fe+1) . 

Compute the QR-factorization of U k , namely, U k = Q k R k , 
Q k £C Nx ( k+1 \ R t ec(' !+1 > x (' ! + 1 ), 

Q k unitary, R k upper triangular with positive diagonal: 


^00 


Q k = [qo\qi I ••• Rk = 

QiQj = $ij V i,j\ r ij = q*u j Vi,j, 


7n • • • r 0ik 
n 1 ••• ri,k 

T'k,k 

ru> 0 V i. 


Step 2. Compute the SVD of R kl namely, R k = YT,H*, 

Y,H G C ( ~ k+1 '> x( ~ k+1 '> unitary, £ £ j£(fc+i)x(fc+i) diagonal: 
Y=[y 0 \Vi\-- - \ y k ], H = [h 0 \h 1 \---\h k ], Y*Y = H*H = J, 
£ = diag (< 7 o,cri,... ,a k ), a 0 > o\ > • • • > a k > 0 . 

Step 3. Set c = [co, ci,..., c k ] T = h k , and compute a = Xo=o c i- 
Set 7 ; = Ci/a, i = 0,1,..., k, provided a^O. 

Compute £= [£o>£i)--*>£k-i] T via 

^o = l-7o; 0 = 0-i _ 7j. j = 1 . 

Step 4. Compute s Tti fc via s nj fc = *„ + Q k _ 1 R k -i$, as follows: 

Compute J 7 = i? fe _i|; T 7 = [?y 0 , r?i,..., ? 7 fe-i] T - 
Compute a n> k = x n + Q A ._ x ry = £c„ + X^do 


Table 4.1: Algorithm for SVD-MPE. 


5 Determinant representations for SVD-MPE 

In [H3 and [23], determinant representations were derived for the vectors s Ujk that 
are produced by the vector extrapolation methods MPE, RRE, MMPE, and TEA. 
These representations have turned out to be very useful in the analysis of the algebraic 
and analytic properties of these methods. In particular, they were used for obtaining 
interesting recursion relations among the s nyk and in proving sharp convergence and 
stability theorems for them. We now derive two analogous determinant representations 
for s U) k produced by SVD-MPE. 

The following lemma, whose proof can be found in [231 Section 3], will be used in 
this derivation in Theorem 15.21 

Lemma 5.1 Let and 7 j be scalars and let the 7 j satisfy the linear system 

k 

^2 u h fi 3 = 0 , i = 0 , 1 ,..., k — 1 , 

T ^ 

= L 

j =0 
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Then, whether Vj are scalars or vectors, there holds 


where 


k 

3=0 


D{y o,V!,...,v k ) = 


? • • • 

D{ 1,1,..., 

1) ’ 



(5.2) 

vo 

Vl 


Vk 



Uo ,0 

uo,i 


U 0 ,k 



u 1,0 

«1,1 


^1 ,k 

5 

(5.3) 

Uk- 1,0 

Uk- 1,1 


Uk-l,k 



Vi are vectors, 

the determinant D{y o, v \,. 

• •, v k ) 


is defined via its expansion with respect to its first row. 


For convenience of notation, we will write 


J-'k 


O 

o 

loi ■ • 

l()k 

^10 

Zn • • 

' Ilk 

IkO 

Ikl ■ ■ 

Ikk 


Then Lk~\ is the principal submatrix of Lf z obtained by deleting the last row and the 
last column of Lk . In addition, Lk-\ is hermitian positive definite just like Lk- 

Theorem 15.21 that follows gives our first determinant representation for s n k result¬ 
ing from SVD-MPE and is based only on the smallest singular value Okk of Uk and 
the corresponding right singular vector hkk- 


Theorem 5.2 Define 

Uij = ( , u n _|_j, u n j r j ')m <Jkk{Lk)iji i,j = 0,1,..., k, (5-4) 


and assume that 

a k k < Ofc-gfc-iH (5.5) 

Then, provided D{ 1,1,... , 1) 0, s nj k exists and has the determinant representation 


*n,k — 


D(x n , ®n+l i • • ■ i ®n+fc) 
D{ 1 , 1 ,...,!) : 


(5.6) 


where D(v o, vi, ..., Vk) is the (fc+1) x (fe+l) determinant defined as in (15.31) in Lemma 
15. II with the Uij as in m - 


Proof. With Uk as in (12.161) . we start by rewriting (|2.18D in the form 


{UkUk — c r kkl)hkk — 0. (5.7) 

3 From the Cauchy interlace theorem, we already know that akk < crk-i,k-i- 
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Invoking here h kk = L k ~c, which follows from (12.191) . and multiplying the resulting 
equality on the left by L f . , we obtain 

(I U* k MU k - a 2 kk L k )c = 0. (5.8) 

Dividing both sides of this equality by j=o c j > an d invoking (12.13ft , we have 

(U* k MU k - a 2 kk L k ) 7 = 0, (5.9) 

which, by the fact that 


(UkMU^j — — (rin-i-j, u n -\.j ) jw, 

is the same as 

k k 

^ '\ { u n+ii u n+j') M ~ (T kk(J- , k)ij\'Yj = 0 ^ 1 u iJ n fj = 0) * = 0, 1, . . . , fc, (5.10) 

i=o i=o 

where we have invoked (15.41) . We will be able to apply Lemma [5. II to prove the validity 
of (15.61) if we show that, in (15.101) . the equations with i = 0,1,..., k — 1, are linearly 
independent, or, equivalently, the first k rows of the matrix 

B k = U* k MU k - a 2 kk L k 


are linearly independent. By the fact that 
Uk \U k — 11 u n _\_ k ] and L k 

we have 


■-1 

hi 

?r 

1 

lk 

/* 

L L k 

Ikk 


j lk [^0fcj Wki • • ■ j lk— l,fc] i 


B k = 


r B ' k -1 

p 

L p * 

p \ 


where 

and 


B’k-i — U* k -iMU k _i — cr kk L k -i 


P — U^MUn+k U kk lki P — u n+k-^- U n+k ^kk^kk- 
Invoking U k _\ = M~ 1/,2 U k -iL 1 J^ l , we obtain 

B'k-i = L^ul^Uk-i - 4kI)Ll^ v 


C\ T* 

Since cr^_ 1 fc _ 1 is the smallest eigenvalue of U k _ 1 U k ^i and since a k -i,k-i > &kk, it 
turns out that B' k _^ is positive definite, which guarantees that the first k rows of B k 
are linearly independent. This completes the proof. ■ 


Remark: We note that the condition that D( 1,1... , 1) ^ 0 in Theorem 15.21 is equiv¬ 
alent to the condition that o c j 7^ 0) which we have already met in Section [2j 


The determinant representation given in Theorem 15.31 that follows is based on the 
complete singular value decomposition of U k , hence is different from that given in 
Theorem 15.21 Since there is no room for confusion, we will denote the singular values 
o k i and right and left singular vectors h ki and g ki of U k by cq, hi and respectively. 
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Theorem 5.3 Let U k be as in (12.161) . and let 


U k = GT,H*, Ge C Nx ( k+1 \ H g c (fc+1)x(fc+1) , S6R (H1)x(t+1) 

be the singular value decomposition ofU k ; that is, 

G = [g 0 \9i | • • • \g k ], 9i9j = $ij‘, H = [h 0 \h 1 \ ■■■ \h k ], h*h j = 5 ij , 


and 

£ = diag (<t 0 , Cl,..., cr k ), a 0 > <ti > ... > a k . 

Define 

u i,j = {M 1/2 9i)*u n + j = 9 *M 1/2 u n+j , i = 0, 1 ,. .. , k - 1 , j = 0,l,...,k, (5.11) 
Then, has the determinant representation 

D(x n , • • • 5 ^n+fc) 




^( 1 , 1 , ■■•,!) 


(5.12) 


where D(v o, ui,..., Ufc) is the (k+ 1) x (k+ 1) determinant defined as in (I5.3j) in Lemma 
5.11 with the Uij as in (15.111) . 


Proof. By Theorem ll.il 

U k h k = a k g k and g*g k = 0, i = 0,1,..., k - 1 . (5.13) 


Therefore, 


g*U k h k = 0, 1 = 0, l,...,fc-l . 


(5.14) 


_-| /q 

By (12.161) and by the fact that c = ' h/-, which follows from (12.191) . and by the fact 

that 7 = c/a, a = ^7=01: which follows from (12.131) . and by (15.141) . we then have 


g\M x l 2 U kl = a~ 1 {g*M 1 ^ 2 U k c) = a~\g*U k h k ) = 0, 1 = 0,1 ,..., A: - 1 . (5.15) 


But, by (15.111) . (15.151) is the same as 

k 

Y. Ui,jlj = 0; 1 = 0,1 ,... ,k - 1. 

1=0 

Therefore, Lemma [5J] applies with Uij as in (15.111) . and the result follows. ■ 


6 SVD-MPE for linearly generated sequences 

In Sidi [l7| , we discussed the connection of the extrapolation methods MPE, RRE, and 
TEA with Krylov subspace methods. We now want to extend the treatment of m to 
SVD-MPE. Here we recall that a Krylov subspace method is also a projection method 
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and that a projection method is defined uniquely by its right and left subspacesQ In 
the next theorem, we show that SVD-MPE is a bone fide Krylov subspace method and 
we identify its right and left subspaces. 

Since there is no room for confusion, we will use the notation of Theorem 15.31 

Theorem 6.1 Let s be the unique solution to the linear system Cx = d, which we 
express in the form 

(J — T)x = d =>- x = Tx + d] T = I — C, 
and let the vector sequence {x m } be produced by the fixed-point iterative scheme 

T d,, Tit — 0,1,... . 

Define the residual vector of x via r{x) = d — Cx. Let also s k = Sg.fc be the approxi¬ 
mation to s produced by SVD-MPE. Then the following are true: 

1. s k is of the form 


k-i 

s k = x 0 + ^dj(CV 0 ) for some dp, r 0 = r(x 0 ) = d - Cx 0 . (6.1) 

i=0 

2. The residual vector of s k , namely, r(s k ), is orthogonal to 
span {M 1/2 g 0 , M l/2 g 1 ,..., M 1/2 g k _ 1 }. Thus, 

(M * 1 / 2 g i )*r(s k ) = 0, * = 0,1,..., k - 1. (6.2) 

Consequently, SVD-MPE is a Krylov subspace method for the linear system Cx = d, 
with the Krylov subspace JC k (C',ro) = spanjro, Ctq, ..., C k ~ 1 ro} as its right subspace 
and span {M L/2 g 0 , M L/2 g 1 ,..., M 1/2 g k _ 1 } as its left subspace. 

Proof. With the x m generated as above, we have 

'U j m +1 — Tu m y — T U o, 771 — 0, 1, . . . . 

Therefore, 

wo = x\ — xq = Tx o + d — xq = d — Cx o = t(xq) 

and 

k —1 k —1 k —1 

&k = X 0 + ^ ^ ~ X 0 + ^ ^ £iT'u o — Xq + ^ ^ Vq. 

2—0 2—0 2—0 

Upon substituting T = I — C in this equality, we obtain (16.11) . 

To prove (16.21) . we first recall that U k ~f = r(s k ) by (13.11) . By this and by (|5.15l) . 
the result in (16.21) follows. ■ 

4 A projection method for the solution of the linear system Cx = d, where C £ C Nx , j g q e g nec j 
as follows: Let y and Z be fc-dimensional subspaces of C N and let xo be a given vector in C N . 
Then the projection method produces an approximation s t, to the solution of Cx = d as follows: 

(i) Sfc = xo + y, y £ T, (h) h*r(sk) = 0 for every h £ Z. y and Z are called, respectively, the right 

and left subspaces of the method. If T is the Krylov subspace )Ck{C; r o) = spanfro, Cro ,..., C k ~ 1 ro}, 
where ro = d — Cxq, then the projection method is called a Krylov subspace method. 
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7 Numerical examples 


We now provide two examples that show the performance of SVD-MPE and compare it 
with MPE. In both examples, SVD-MPE and MPE are implemented with the standard 
Euclidean inner product and the norm induced by it. Thus, M = I and Lk = I 
throughout. 

As we have already mentioned, a major application area of vector extrapolation 
methods is that of numerical solution of large systems of linear or nonlinear equations 
?/>( x) = 0 by fixed-point iterations x m+ i = f(x m ). [Here x = f(x) is a possibly 
preconditioned form of i/)(x) = 0.] For SVD-MPE, as well as all other polynomial 
methods discussed in the literature, the computation of the approximation s n ^ to s, 
the solution of i/)(x) = 0, requires k +1 of the vectors x m to be stored in the computer 
memory. For systems of very large dimension N, this means that we should keep A: at a 
moderate size. In view of this limitation, a practical strategy for systems of equations 
is cycling, for which n and k are fixed. Here are the steps of cycling: 

CO. Choose integers n > 0, and k > 1, and an initial vector xq. 

Cl. Compute the vectors X\,X 2 , ■ ■ ■ ,x n+ k+i [via x m+ i = f(x m )]. 

C2. Apply SVD-MPE (or MPE) to the vectors x n , x n+ i ,..., x n+ k+\, with end result 

Sn,k ■ 

C3. If s nt k satisfies accuracy test, stop. 

Otherwise, set xq = s n and go to Step Cl. 

(r) 

We will call each application of steps C1-C3 a cycle, and denote by s n ), the s n ^ that is 
computed in the rth cycle. We will also denote the initial vector xq in step CO by s^\,. 

Numerical examples suggest that the sequence {s^,}£T 0 has very good convergence 
properties. 


Example 7.1 Consider the vector sequence {x m } obtained from x m+ i = Tx m + d, 
m = 0 , 1 ,..., where 


T = 0.06 x 


5 2 1 

2 6 3 

1 3 6 

1 1 3 

1 1 


1 

1 1 
3 1 

6 3 


and is symmetric with respect to both main diagonals, and T € H NxN The vector d 
is such that the exact solution to x = Tx + d is s = [1,1,... , 1] T . We have p(T) < 1, 
so that {x m } converges to s. 

Figure [T] shows the I 2 norms of the errors in s n *., n = 0,1,..., with k = 5 fixed. 
Here N = 100. Note that all of the approximations s n< 5 make use of the same (infinite) 
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vector sequence { x m }, and, practically speaking, we are looking at how the methods 
behave as n —>• oo. It is interesting to see that SVD-MPE and MPE behave almost 
the same. Although we have a rigorous asymptotic theory confirming the behavior of 
MPE in this mode as observed in Figured] (see }16] . m-m), we do not have any 
such theory for SVD-MPE at the present. 

Figure [ 2 ] shows the I 2 norm of the error in s n & in the cycling mode with n = 0 and 
k = 20. Now N = 1000. 

Example 7.2 We now apply SVD-MPE and MPE to the nonlinear system that arises 
from finite-difference approximation of the two-dimensional convection-diffusion equa¬ 
tion 


—V 2 u + Cu(u x + u y ) = /, (x, y) € D = ( 0 ,1) x ( 0 ,1), 

where u(x,y ) satisfies homogeneous boundary conditions. f(x,y) is constructed by 
setting C = 20 in the differential equation and by taking 

u(x, y) = 10 xy(l — x)(l — y) exp(x 4 " 5 ) 

as the exact solution. 

The equation is discretized on a square grid by approximating u xx , u yy , u x , and u y 
by centered differences with truncation errors 0(h 2 ). Thus, letting h = \/v, and 

(xuyj) = (■ ih,jh ), i,j = 0, l,...,zz, 


and 


u x (xi,yj ) 


u(x i+1 ,yj) - u(xj -1 ,yj) 
2 h 


", Uy 




u(xj,y j+1 ) - u{xj,yj- 1 ) 

2 h 


and 


- V 2 u(xi,yj) 


4 u(xi,yj) - u(x i+ i,yj) - u(xi-i,yj) - u(xi , y j+1 ) - u(xi,yj- 1 ) 

h 2 


we replace the differential equation by the finite difference equations 


4 Uij Ui -|_ij Ui—ij 1 Ujj—\ 

h 2 

1 ri„. u i+l,j ~ u i~l,j + u i,j+l ~ u i,j~ 1 
+ Cmj 2h 


f(xi,yj), l<i,j <v-l, 


with 

UO j — ^i ,0 — — 'Wj,!/ — 0 V i, j. 

Here is the approximation to u(xi,yj), as usual. 

We first write the finite difference equations in a way that is analogous to the PDE 
written in the form 

—V 2 it = / - Cu(u x + Uy), 
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and split the matrix representing —V 2 to enable the use of the Jacobi and Gauss-Seidel 
methods as the iterative procedures to generate the sequences {x m }. 

Figures [3] and |4] show the I 2 norms of the errors in s n j ; , from SVD-MPE and MPE in 
the cycling mode with n = 0 and k = 20, the iterative procedures being, respectively, 
that of Jacobi and that of Gauss-Seidel for the linear part —V 2 u of the PDE. Here 
u = 100, so that the number of unknowns (the dimension) is N = 99 2 . 


Acknowledgement 

The author would like to thank Boaz Ophir for carrying out the computations reported 
in Section [7] of this work. 


References 

[1] C. Brezinski. Application de 1’e-algorithme a la resolution des systenres non 
lineaires. C. R. Acad. Sci. Paris , 271 A:1174-1177, 1970. 

[2] C. Brezinski. Sur un algorithme de resolution des systemes non lineaires. C. R. 
Acad. Sci. Paris , 272 A:145-148, 1971. 

[3] C. Brezinski. Generalisations de la transformation de Shanks, de la table de Pade, 
et de l’e-algorithme. Calcolo , 11:317-360, 1975. 

[4] S. Cabay and L.W. Jackson. A polynomial extrapolation method for finding limits 
and antilimits of vector sequences. SIAM J. Numer. Anal, 13:734-752, 1976. 

[5] R.P. Eddy. Extrapolating to the limit of a vector sequence. In P.C.C. Wang, 
editor, Information Linkage Between Applied Mathematics and Industry, pages 
387-396, New York, 1979. Academic Press. 

[6] E. Gekeler. On the solution of systems of equations by the epsilon algorithm of 
Wynn. Math. Comp., 26:427-436, 1972. 

[7] G.H. Golub and W. Kahan. Calculating the singular values and pseudo-inverse 
of a matrix. SIAM J. Numer. Anal, Series B, 2:205-224, 1965. 

[8] G.H. Golub and C.F. Van Loan. Matrix Computations. The Johns Hopkins 
University Press, Baltimore, fourth edition, 2013. 

[9] P.R. Graves-Morris. Vector valued rational interpolants I. Numer. Math., 42:331- 
348, 1983. 

[10] P.R. Graves-Morris. Extrapolation methods for vector sequences. Numer. Math., 
61:475-487, 1992. 

[11] R.A. Horn and C.R. Johnson. Matrix Analysis. Cambridge University Press, 
Cambridge, 1985. 


18 


[12] A.S. Householder. The Theory of Matrices in Numerical Analysis. Blaisedell, New 
York, 1964. 

[13] M. Mesina. Convergence acceleration for the iterative solution of the equations 
X = AX +/. Comput. Methods Appl. Mech. Engrg ., 10:165-173, 1977. 

[14] B.P. Pugachev. Acceleration of the convergence of iterative processes and a 
method of solving systems of nonlinear equations. U.S.S.R. Comput. Math. Math. 
Phys., 17:199-207, 1978. 

[15] D. Shanks. Nonlinear transformations of divergent and slowly convergent se¬ 
quences. J. Math, and Phys., 34:1-42, 1955. 

[16] A. Sidi. Convergence and stability properties of minimal polynomial and reduced 
rank extrapolation algorithms. SIAM J. Numer. Anal., 23:197-209, 1986. Origi¬ 
nally appeared as NASA TM-83443 (1983). 

[17] A. Sidi. Extrapolation vs. projection methods for linear systems of equations. J. 
Comp. Appl. Math., 22:71-88, 1988. 

[18] A. Sidi. Efficient implementation of minimal polynomial and reduced rank extrap¬ 
olation methods. J. Comp. Appl. Math., 36:305-337, 1991. Originally appeared 
as NASA TM-103240 ICOMP-90-20 (1990). 

[19] A. Sidi. Convergence of intermediate rows of minimal polynomial and reduced 
rank extrapolation tables. Numer. Algorithms, 6:229-244, 1994. 

[20] A. Sidi. Vector extrapolation methods with applications to solution of large sys¬ 
tems of equations and to PageRank computations. Comp. & Maths, with Applies., 
56:1-24, 2008. 

[21] A. Sidi. Review of two vector extrapolation methods of polynomial type with 
applications to large-scale problems. J. Comput. Sci., 3:92-101, 2012. 

[22] A. Sidi and J. Bridger. Convergence and stability analyses for some vector extrap¬ 
olation methods in the presence of defective iteration matrices. J. Comp. Appl. 
Math., 22:35-61, 1988. 

[23] A. Sidi, W.F. Ford, and D.A. Smith. Acceleration of convergence of vector se¬ 
quences. SIAM J. Numer. Anal., 23:178-196, 1986. Originally appeared as NASA 
TP-2193 (1983). 

[24] A. Sidi and Y. Shapira. Upper bounds for convergence rates of acceleration meth¬ 
ods with initial iterations. Numer. Algorithms, 18:113-132, 1998. 

[25] D.A. Smith, W.F. Ford, and A. Sidi. Extrapolation methods for vector sequences. 
SIAM Rev., 29:199-233, 1987. Erratum: SIAM Rev., 30:623-624, 1988. 

[26] J. Stoer and R. Bulirsch. Introduction to Numerical Analysis. Springer-Verlag, 
New York, third edition, 2002. 


19 



[27] L.N. Trefethen and D. Bau, III. Numerical Linear Algebra. SIAM, Philadelphia, 
1997. 

[28] P. Wynn. On a device for computing the e m (S n ) transformation. Mathematical 
Tables and Other Aids to Computation , 10:91-96, 1956. 

[29] P. Wynn. Acceleration techniques for iterated vector and matrix problems. Math. 
Comp., 16:301-322, 1962. 

[30] P. Wynn. Continued fractions whose coefficients obey a noncommutative law of 
multiplication. Arch. Rat. Mech. Anal., 12:273-312, 1963. 

[31] P. Wynn. General purpose vector epsilon algorithm procedures. Numer. Math., 
6:22-36, 1964. 


20 



Error Norm 



Figure 1: I 2 norm of error in s n k, n = 0,1,, with k = 5, from MPE and SVD-MPE, 
for Example 17.11 with N = 100. 
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Error Norm 



Figure 2: I 2 norm of error in So,20 i n the cycling mode, from MPE and SVD-MPE, for 
Example 17.II with N = 1000. 
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Jacobi 



Figure 3: I 2 norm of error in So,20 i n the cycling mode, from MPE and SVD-MPE, for 
Example 17.21 with u = 100 hence N = 99 2 . The underlying iteration method is that of 
Jacobi. 
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Gauss-Seidel 



Number of Iterations 

Figure 4: I 2 norm of error in So,20 i n the cycling mode, from MPE and SVD-MPE, for 
Example 17.21 with u = 100 hence N = 99 2 . The underlying iteration method is that of 
Gauss-Seidel. 
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